{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kriging in python\n",
    "\n",
    "Here, ordinary kriging using the geostat package is illustrated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import geostat as geo\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A data set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAEk9JREFUeJzt3XGMpPdd3/H3x1wT9UyCU24DcXyX\nDeh8SKTOYRaa/kFxYiU4TnWhMhGkJ4xcqsUQoxA1TYxWtStZpyaYyA0KCTolhwGtXYXUhWBDsEkb\nTEVNu5falwtxHTfcnQ8nuTVubSmnEjn+9o95zp5bdnf2ZmdvZn73fkmrmec7z+x+tNr53HPPzPwm\nVYUkqV0XjTuAJGlrWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxm0bdwCAHTt2\n1Ozs7LhjSNJUOXz48FNVNTNov4ko+tnZWZaWlsYdQ5KmSpLjG9nPUzeS1DiLXpIaZ9FLUuMseklq\nnEUvSY2z6CU15Y4HHht3hIlj0Utqyoc/++VxR5g4Fr0kNW4i3jAlSZtxxwOPnXUkP3vzfQC8++rd\nvOfNl48r1sTIJHw4+NzcXPnOWEmjMHvzfRz7wNvGHeO8SHK4quYG7eepG0lqnEUvqSnvvnr3uCNM\nHIteUlM8J/93WfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDVuYNEnOZTkVJKjq9z23iSV\nZEe3nSS/luTxJEeSXLkVoafZH978IZidhYsu6l0uLo47kjSQS/9Ot40c0d8JXLNymGQn8GbgRN/4\nrcDu7mse+NjmIzZkcZGrPrQAx49DVe9yft6y18Rz6d/pNrDoq+pB4OlVbroDeB/Qvyra24Hfrp6H\ngEuSvGokSVuwsMD25/727Nnp07CwMJ48ki4IQy1TnGQf8NdV9UiS/pteDTzRt32ym311le8xT++o\nn127dg0TY2qcWUL1K8dPrP4v64kTq02lsXLp33ZsaJniJLPAvVX1uiTbgf8CvKWqnklyDJirqqeS\n3Af8u6r6r939Pgu8r6oOr/f9L5hlimdne6drVnrNa+DYsfOdRtqwC2np32mylcsUfy/wWuCRruQv\nAz6f5LvpHcHv7Nv3MuDJIX5Gmw4c4PS2l549274dDhwYTx5JF4RzLvqq+kJVvbKqZqtqll65X1lV\nXwM+DVzfvfrmDcAzVfV3TttcsPbv53P/6kDvCD7pXR48CPv3jzuZtC6X/p1uA0/dJLkbuArYAXwd\nuLWqPtF3+zFePHUT4CP0XqVzGrihqgaek7lgTt1I0ght9NTNwCdjq+qdA26f7btewLs2ElCSdH74\nzlhJapxFL0mNs+glqXEWvSQ1zqKXpMZZ9JLUOItekho39UXvOtmStL6pL3rXyZak9U190UuS1jfU\nevTj5jrZkrRxG1qPfqttZlEz18mWdKHayvXoJUlTZOqL3nWyJWl9U1/0npOXpPVNfdFLktZn0UtS\n4yx6SWqcRS9JjbPoJalxFr0kNW5g0Sc5lORUkqN9s9uSHEnycJL7k1zazb8jyR8keSTJF5PcsJXh\nJUmDbeSI/k7gmhWz26vqiqraC9wL3NLN3wX8ZVW9HrgK+FCSl4woqyRpCAOLvqoeBJ5eMXu2b/Ni\n4MyCOQW8LEmAb+/u99xookqShjH06pVJDgDXA88Ab+zGHwE+DTwJvAz4yap6frMhJUnDG/rJ2Kpa\nqKqdwCJwUzf+MeBh4FJgL/CRJC9f7f5J5pMsJVlaXl4eNoYkaYBRvOrmLuC67voNwD3V8zjwV8D3\nrXanqjpYVXNVNTczMzOCGJKk1QxV9En6l4zcBzzaXT8BXN3t813AHuArmwkoSdqcgefok9xN7xU0\nO5KcBG4Frk2yB3geOA7c2O1+G3Bnki8AAd5fVU9tRXBJ0sYMLPqqeucq40+sse+TwFs2G0qSNDq+\nM1aSGmfRS1LjLHpJapxFL0mNs+glqXEWvSQ1zqKXpMZZ9JLUOItekhpn0UtS4yx6SWqcRS9JjbPo\nJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDVuYNEnOZTkVJKjfbPb\nkhxJ8nCS+5Nc2nfbVd38i0n+dKuCS5I2ZiNH9HcC16yY3V5VV1TVXuBe4BaAJJcAHwX2VdX3A+8Y\nYVZJ0hAGFn1VPQg8vWL2bN/mxUB11/85cE9Vnej2OzWinJKkIQ19jj7JgSRPAPvpjuiBy4FXJPlc\nksNJrl/n/vNJlpIsLS8vDxtDkjTA0EVfVQtVtRNYBG7qxtuAHwTeBvwY8G+SXL7G/Q9W1VxVzc3M\nzAwbQ5I0wChedXMXcF13/STwmar6RlU9BTwIvH4EP0OSNKShij7J7r7NfcCj3fXfB34kybYk24F/\nBHxpcxElSZuxbdAOSe4GrgJ2JDkJ3Apcm2QP8DxwHLgRoKq+lOQzwJHuto9X1dFVv7Ek6bxIVQ3e\na4vNzc3V0tLSuGNI0lRJcriq5gbt5ztjJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuMseklqnEUv\nSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXEWvSQ1zqKXpMZZ9JLU\nOItekho3sOiTHEpyKsnRvtltSY4keTjJ/UkuXXGfH0ryrSQ/sRWhJUkbt5Ej+juBa1bMbq+qK6pq\nL3AvcMuZG5J8G/BB4I9HFVJqxuIizM7CRRf1LhcXx51IF4CBRV9VDwJPr5g927d5MVB9278I/Efg\n1CgCSs1YXIT5eTh+HKp6l/Pzlr223NDn6JMcSPIEsJ/uiD7Jq4F/BvzGaOJJDVlYgNOnz56dPt2b\nS1to6KKvqoWq2gksAjd1438PvL+qvjXo/knmkywlWVpeXh42hjQ9Tpw4t7k0IqN41c1dwHXd9Tng\nPyQ5BvwE8NEkP77anarqYFXNVdXczMzMCGJIE27XrnObSyMyVNEn2d23uQ94FKCqXltVs1U1C3wK\n+IWq+r1Np5RacOAAbN9+9mz79t5c2kLbBu2Q5G7gKmBHkpPArcC1SfYAzwPHgRu3MqTUhP37e5cL\nC73TNbt29Ur+zFzaIqmqwXttsbm5uVpaWhp3DEmaKkkOV9XcoP18Z6wkNc6il6TGWfSS1DiLXpIa\nZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXEW\nvSQ1zqKXpMZZ9JLUOItekhpn0UtS4wYWfZJDSU4lOdo3uy3JkSQPJ7k/yaXdfH83P5Lkz5O8fivD\nS5IG28gR/Z3ANStmt1fVFVW1F7gXuKWb/xXwo1V1BXAbcHBUQdWAxUWYnYWLLupdLi6OO5F0Qdg2\naIeqejDJ7IrZs32bFwPVzf+8b/4QcNnmI6oJi4swPw+nT/e2jx/vbQPs3z++XNIFYOhz9EkOJHkC\n2M+LR/T9fhb4o2G/vxqzsPBiyZ9x+nRvLmlLDV30VbVQVTuBReCm/tuSvJFe0b9/rfsnmU+ylGRp\neXl52BiaFidOnNtc0siM4lU3dwHXndlIcgXwceDtVfU3a92pqg5W1VxVzc3MzIwghibarl3nNpc0\nMkMVfZLdfZv7gEe7+S7gHuCnq+qxzcdTMw4cgO3bz55t396bS9pSA5+MTXI3cBWwI8lJ4Fbg2iR7\ngOeB48CN3e63AN8JfDQJwHNVNbcFuTVtzjzhurDQO12za1ev5H0iVtpyqapxZ2Bubq6WlpbGHUOS\npkqSwxs5mPadsZLUOItekhpn0UtS4yx6SWqcRS9JjbPoJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FL\nUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXEDiz7J\noSSnkhztm92W5EiSh5Pcn+TSbp4kv5bk8e72K7cyvCRpsI0c0d8JXLNidntVXVFVe4F7gVu6+VuB\n3d3XPPCxEeWUJA1pYNFX1YPA0ytmz/ZtXgxUd/3twG9Xz0PAJUleNaqwkqRzt23YOyY5AFwPPAO8\nsRu/Gniib7eT3eyrw/4cSdLmDP1kbFUtVNVOYBG4qRtntV1Xu3+S+SRLSZaWl5eHjSFJGmAUr7q5\nC7iuu34S2Nl322XAk6vdqaoOVtVcVc3NzMyMIIYkaTVDFX2S3X2b+4BHu+ufBq7vXn3zBuCZqvK0\njSSN0cBz9EnuBq4CdiQ5CdwKXJtkD/A8cBy4sdv9D4FrgceB08ANW5BZknQOBhZ9Vb1zlfEn1ti3\ngHdtNpQkaXR8Z6wkNc6il6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6BtxxwOPjTuCpAll0Tfi\nw5/98rgjSJpQFr0kNW7o9eg1fnc88NhZR/KzN98HwLuv3s173nz5uGJJmjDpLU8zXnNzc7W0tDTu\nGFNt9ub7OPaBt407hqTzKMnhqpobtJ+nbiSpcRZ9I9599e7BO0m6IFn0jfCcvKS1WPSS1DiLXpIa\nZ9FLUuMseklqnEUvSY2z6CWpcQOLPsmhJKeSHO2b3Z7k0SRHkvynJJd087+X5LeSfCHJl5L88laG\nlyQNtpEj+juBa1bMHgBeV1VXAI8BZwr9HcBLq+ofAj8I/FyS2ZEklcbEJaC1lc7H39fAoq+qB4Gn\nV8zur6rnus2HgMvO3ARcnGQb8PeBbwLPji6udP65BLS20vn4+xrFOfp/AfxRd/1TwDeArwIngF+t\nqqfXuqMkaettapniJAvAc8BiN/ph4FvApcArgD9L8idV9ZVV7jsPzAPs2rVrMzGkkXMJaG2l8/33\ntaFlirvz7PdW1ev6Zj8D3AhcXVWnu9mvAw9V1e9024eAz1TVJ9f7/i5TrEnmEtDaSpv5+9rSZYqT\nXAO8H9h3puQ7J4A3pedi4A3Ao8P8DEnSaGzk5ZV3A/8N2JPkZJKfBT4CvAx4IMnDSX6j2/3XgW8H\njgL/A/jNqjqyNdGl88MloLWVzsffl58wJUlTyk+YkiQBFr0kNc+il6TGWfSS1DiLXpIaNxGvukmy\nDBwfd44VdgBPjTvEOZimvNOUFaYr7zRlhenKO4lZX1NVM4N2moiin0RJljbysqVJMU15pykrTFfe\nacoK05V3mrKu5KkbSWqcRS9JjbPo13Zw3AHO0TTlnaasMF15pykrTFfeacp6Fs/RS1LjPKKXpMZZ\n9J0klyT5VPeh519K8o/X+hD0cVsta99t701SSXaMM+MZa2VN8otJ/leSLyb5lXHnPGONv4O9SR7q\nVmpdSvLD484JkGRPl+nM17NJfinJP0jyQJIvd5evmOCsk/oYWzVv3+0T9TgbqKr86p2++i3gX3bX\nXwJcArwF2NbNPgh8cNw518raXd8J/DG99yTsGHfOdX6vbwT+hN4HyQO8ctw5B+S9H3hrN7sW+Ny4\nc66S+9uArwGvAX4FuLmb3zwpf7drZJ3Ix9haebvtiXucDfryiB5I8nLgnwCfAKiqb1bV/621PwR9\nbNbK2t18B/A+eh/SPnbrZP154ANV9bfd/NT4Ur5onbwFvLzb7TuAJ8eTcF1XA/+7qo4Db6f3Dxbd\n5Y+PLdXqXsg6iY+xVfT/bmHCHmcbYdH3fA+wDPxmkv+Z5OPdJ2T16/8Q9HFaNWuSfcBfV9UjY87X\nb63f6+XAjyT5iyR/muSHxhvzBWvl/SXg9iRPAL8K/PI4Q67hp4C7u+vfVVVfBeguXzm2VKvrz9pv\nUh5jK72Qd0IfZwNZ9D3bgCuBj1XVDwDfoPdfXmDVD0Efp9Wy/ltgAbhljLlWs9bvdRu9D49/A/Cv\ngU8mydhSvmitvD8PvKeqdgLvoTvinxRJXgLsA3533FkGWSvrhD3GXtCfN8l2JvNxNpBF33MSOFlV\nf9Ftf4reA/7Mh6D/U2B/dSfoxmytrK8FHklyjN5/fz+f5LvHE/EFa2U9CdxTPf8deJ7eOiLjtlbe\nnwHu6Wa/C0zEk7F93gp8vqq+3m1/PcmrALrLiTg11lmZdRIfY/36834vk/k4G8iiB6rqa8ATSfZ0\no6uBv1znQ9DHZo2sn6+qV1bVbFXN0iusK7t9x2at3yvwe8CbAJJcTu9Jz7EvFrVO3ieBH+1mbwK+\nPIZ463knZ58K+TS9f5zoLn//vCda21lZJ/ExtsILeavqC5P4ONsI3zDVSbIX+Di90vkKcAO9Dzh/\nKfA33W4PVdWN40n4otWyVtX/6bv9GDBXVWMvzzV+r98ADgF7gW8C762q/zy2kH3WyPv9wIfpndr5\nf8AvVNXhsYXs051OeAL4nqp6ppt9J/BJYBdwAnhHVT09vpQ9a2R9nAl8jMHqeVfcfowJeZwNYtFL\nUuM8dSNJjbPoJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuMseklq3P8Hkkqc/KkFJs8AAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1dfef0f3518>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Input data: the measurements\n",
    "x1 = np.array( [61,63,64,68,71,73,75] )        # Locations of the points (x1=x, x2=y)\n",
    "x2 = np.array( [139,140,129,128,140,141,128] )  \n",
    "v = np.array( [477,96,227,646,606,791,783] )   # Value of the measurements \n",
    "\n",
    "# The data in matrix form = one line per data point\n",
    "x = np.transpose( (x1,x2) )\n",
    "\n",
    "# The locations where we want to estimate \n",
    "x10 = np.array([63.1, 66, 69]); \n",
    "x20 = np.array([140, 132, 134]);\n",
    "xu = np.transpose( (x10,x20) )\n",
    "\n",
    "# Plot the data for visual check\n",
    "plt.plot(x1,x2,'+');   # Position des points\n",
    "plt.plot(x10,x20,'or') # Point où l'on veut estimer\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kriging using the geostat module\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 102.36352305]\n",
      " [ 395.41363864]\n",
      " [ 526.231114  ]]\n",
      "[[ 39.78788501]\n",
      " [ 33.88010645]\n",
      " [ 32.17473703]]\n"
     ]
    }
   ],
   "source": [
    "# Defines the type of covariance function and parameters\n",
    "import geostat as geo\n",
    "covar = geo.covariance(rang=30, sill=20, typ=\"spherical\")\n",
    "vo,so = geo.ordinary_kriging_covariance(x,v,xu,covar)\n",
    "print(vo)\n",
    "print(so)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 102.70760736]\n",
      " [ 400.38896026]\n",
      " [ 529.99471018]]\n",
      "[[   3.94873012]\n",
      " [  87.85744942]\n",
      " [ 111.39842006]]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "vario = geo.variogram(rang=30, sill=20, typ=\"spherical\")\n",
    "vo,so = geo.ordinary_kriging_variogram(x,v,xu,vario)\n",
    "print(vo)\n",
    "print(so)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  9.7,  47.6,  18.8],\n",
       "       [ 43.8,  24.6,  67.9]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xu.shape=(2,1)\n",
    "np.concatenate((x,xu),axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
